Gas-liquid critical parameters of asymmetric models of ionic fluids 

O.V. Patsahan 1 and T.M. Patsahan 1 

1 Institute for Condensed Matter Physics of the National Academy of Sciences of Ukraine, 

1 Svientsitskii Str., 79011 Lviv, Ukraine 
(Dated: July 18, 2011) 

Abstract 

The effects of size and charge asymmetry on the gas-liquid critical parameters of a primitive 
model (PM) of ionic fluids are studied within the framework of the statistical field theory based 
on the collective variables method. Recently, this approach has enabled us to obtain the correct 
trends of the both critical parameters of the equisize charge-asymmetric PM without assuming 
ionic association. In this paper we focus on the general case of an asymmetric PM characterized by 
the two parameters: hard-sphere diameter-, A = <r+/o"_ and charge, z = q+/\q_\, ratios of the two 
ionic species. We derive an explicit expression for the chemical potential conjugate to the order 
parameter which includes the effects of correlations up to the third order. Based on this expression 
we consider the three versions of PM: a monovalent size-asymmetric PM (A / 1, z = 1), an equisize 
charge-asymmetric PM (A = 1, z / 1) and a size- and charge-asymmetric PM (A / 1, z = 2). 
Similar to simulations, our theory predicts that the critical temperature and the critical density 
decrease with the increase of size asymmetry. Regarding the effects of charge asymmetry, we obtain 
the correct trend of the critical temperature with z, while the trend of the critical density obtained 
in this approximation is inconsistent with simulations, as well as with our previous results found in 
the higher-order approximation. We expect that the consideration of the higher-order correlations 
will lead to the correct trend of the critical density with charge asymmetry. 
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I. INTRODUCTION 



It is well-known, that electrostatic forces determine the properties of various systems: 
physical as well as chemical or biological. In particular, the Coulomb interactions are of 
great importance when dealing with ionic fluids. Ionic fluids include molten salts, elec- 
trolyte solutions and ionic liquids. In most cases the Coulomb interaction is the dominant 
interaction and due to its long-range character it can substantially affect the critical prop- 
erties and the phase behavior of ionic systems. Thus, the investigations dealing with these 
issues are of fundamental interest and of practical importance. 

Over the last fifteen years, both phase diagrams and the critical behavior of the sys- 
tems with dominant Coulomb interactions have been intensively studied by simulations and 
theoretical methods. A theoretical model that demonstrates the phase separation driven 
exclusively by Coulomb forces is a primitive model (PM). In this model, the ionic fluid is 
described as an electroneutral mixture of charged hard spheres immersed in a structure- 
less dielectric continuum. The simplest version of the two-component PM, its symmetrical 
version, is called a restricted primitive model (RPM). A gas-liquid phase transition of the 
RPM is well established. However, over the years the figures for the critical parameters have 
changed substantially. Now there is a good agreement between the recent simulations per- 
formed by different teams. The estimations turn out to be near T* ~ 0.049, p* ~ 0.06 — 0.08 
when the temperature T* and the density p* are in standard dimensionless units [1]. 

As concerns the theory, several theoretical methods have been proposed in which the ion 
association is explicitly taken into account. The main of them are the generalized Debay- 
Hiickel (GDH) theory and the associated mean spherical approximation 2J. These theories 
are based on the addition of the chemical association model of Bjerrum [3j or Ebeling and 
Grigo 4j. The GDH theory (solvated ion-cluster theory with hard-core term) yields the 
following estimations for the critical parameters of the RPM: T* = 0.0557, p* = 0.0261 {jj. 

More recently, the study of the phase behavior of size- and charge-asymmetric PMs 
has been started. The key findings from simulation studies of asymmetric models are as 
follows: the suitably normalized critical temperatures decrease with size and charge asym- 
metry while the critical densities increase with charge asymmetry but decrease with size 



asymmetry 
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. Comparison of simulated critical parameters and 



theoretical predictions for asymmetric models has revealed that several established theories, 
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such as the mean spherical approximation (MSA) and the original DH theory are not ca- 
pable of predicting the trends observed in simulations Q|- Moreover, both the original 
DH theory and the MSA predict no dependence on charge asymmetry in the equisize case. 
The exception are the theories mentioned above that include the association effects explic- 



itly 



except 



171 ]. The trends found from the GDH theory for the critical parameters of an 



equisize (z:l) charge-asymmetric PM as a function of charge asymmetry qualitatively agree 
with simulation data 5|. As regards the size asymmetry, the extensions of the DH theory 
for monovalent size- asymmetric PMs that describe the charge-unbalanced "border zones" 
surrounding each ion lead to the trends of the both critical parameters that qualitatively 
agree with simulation predictions 18]. However, this is true only for modest size asymme- 
tries. More recently, the study of the effects of size and charge asymmetry on the gas-liquid 
phase separation has been started within the field-theoretical description in [l9j. It is found 
that only some of the effects of the size and charge asymmetry are correctly predicted at 
the mean-field (MF) level of the theory. For example, the trend of T* with size asymmetry 
obtained in this approximation is inconsistent with the predictions of simulations. Summa- 
rizing, we can state that even the qualitative theoretical understanding of the issue is not 
quite appropriate. 

In this paper we focus on the issue of the effects of charge and size asymmetry on the gas- 
liquid critical parameters of PMs. To this end, we use the statistical field theory based on 
the method of collective variables (CVs) (see 20] and the references herein). The approach 
allows us to derive the exact functional representation of the grand partition function and 



formulate, on this basis, the perturbation theory, 
known theories were established recently 
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inks between this approach and the other 



recently 



231 ] . We also use the method proposed 



211 ] for the study of the gas-liquid phase diagram of equisize charge-asymmetric 



PMs. The method is based on determining the chemical potential conjugate to the order 
parameter and allows one to take into account the effects of higher-order correlations. Its 
application to an equisize (z:l) charge-asymmetric PM enabled us to obtain the qualitative 
agreement of the trends of T*(z) and p*(z) with simulation findings 2l[. The theory also 
yiel ds the best theoretical quantitative estimates for the critical parameters of the RPM 
[24( . Here we study the general case of size- and charge-asymmetric PMs. 

The layout of the paper is as follows. In Section 2, starting with the Hamiltonian of a 
size- and charge-asymmetric PM we sketch out the main points of the CVs based theory. We 
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analyze the Gaussian approximation of the functional of the grand partition function and 
determine the CV connected with the order parameter for the gas-liquid phase separation. 
In Section 3 we study the gas-liquid critical parameters of asymmetric PMs taking into 
account the correlation effects of higher order. We consider the three versions of PM: an 
equisize PM with charge asymmetry; a monovalent PM with size asymmetry; a size- and 
(2:1) charge-asymmetric PM. We conclude in Section 4. 



II. COLLECTIVE VARIABLES BASED THEORY FOR ASYMMETRIC PMS 



A. Model 



We consider a classical two-component system consisting of N + hard spheres of diameter 
a + carrying a charge q + = zq and N_ hard spheres of diameter cr_ carrying a charge 
g_ = —q. The ions are immersed in a structureless dielectric medium. The system is 
electrically neutral: Yl,a=+ - QaPa = 0, p a = N a /V is the number density of the ath species. 

The pair interaction potential is assumed to be of the following form: 

M^^tfW+ftw, (i) 

where <Papi r ) ls the interaction potential between the two additive hard spheres of diame- 
ters a a and a p. We call the two-component hard sphere system a reference system (RS). 
Thermodynamic and structural properties of RS are assumed to be known. 0^g(V) is the 
Coulomb potential: <^p{r) = q a q$$ c if) , where <p c {r) = l/(Dr), D is the dielectric constant 
and hereafter we put D = 1. The model is characterized by the parameters of size and 
charge asymmetry: 



The fluid is at equilibrium in the grand canonical ensemble. The grand partition function 
(GPF) of the model (CQ) can be written as follows: 



^ l?/J 2^ 11 — /( dr )exp 

Af + >0iV_>0 «=+,- a ' 



a(3 ij 

where the following notations are used: v a is the dimensionless chemical potential, v a = 
(3 ii a — 3 In A a , fi a is the chemical potential of the ath species, (3 is the reciprocal temperature, 



A^ 1 = {2nm a f3~ l /h 2 ) 1 / 2 is the inverse de Broglie thermal wavelength; (dr) is the element 
of configurational space of the particles. 

It is worth noting that the regularization of the potential 0^(V) inside the hard core is 
arbitrary to some extent. For example, different regularizations for the Coulomb potential 
were considered in 
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261 ] . Within the framework of the Gaussian approximation of GPF 



ihe best estimation for the critical temperature is achieved for the optimized regularization 



27l | that leads to the ORPA (MSA). However^this approximation does not work properly in 



the higher orders of the perturbation theory 



(WCA) regularization scheme for 



261 ] . Here we use the Weeks- Chandler- Andersen 



q a qf3/v a i3, r < O afi 
Qaqp/r, r > a af3 . 



(2) 



B. Functional of the GPF of an asymmetric PM. The method of collective vari- 
ables 



Now we use the CVs based theory, developed in 20|] for a multicomponent continuous 
system with short- and long-range interactions in the grand canonical ensemble. As a result, 
the exact functional representation of the GPF for the PM with the interaction potential 
([!]) can be written in the form: 



(dp) (dw) exp (-H[v a ; p a , u a ]) 



(3) 



where the action TC reads as 



P 



a, (3 k a k 



VjJr 



(4) 



Here p^ a = p k a — ip k a is the CV which describes the value of the k-th fluctuation mode of 
the number density of the ath species, the indices c and s denote real and imaginary parts 
of pk,a- ^k,a is conjugate to the CV p^ a and each of p^ a {oJ^ a ) takes all the real values 
from — oo to +oo. (dp) and (duo) are volume elements of the CV phase space 

(dp) = dp , a Yl d Pk,a d Pk, a > ( dw ) = II dw o. a II dw k,« d ^k,a 

a k^O a k^O 

and the product over k is performed in the upper semi-space (p~k,« = p^ a , ^-k,a = k'ka)- 



ap(k) is the Fourier transform of the Coulomb potential (r). In the case of the WCA 

r 0( 



regularization (see ([2])) we obtain for 



T*z(l —5) x 

m _ (k) = -^2iHM, (7) 



where the following notations are introduced: 



T = 8 

is the dimensionless temperature, 

and x = ka±, a± = (cr + + cr_)/2. 

Shs [^a — i^a] is the GPF of a two-component hard sphere system with the renormalized 
chemical potential 

in the presence of the local field — iu a (r). 

In order to formulate the perturbation theory we present the CVs in the following form: 

where the mean-field (MF) values p a and uj a are the solutions of the saddle point equations. 
Then we present lnSnsf^a — iu) a ] in (j3J) in the form of the cumulant expansion 
(— i) n 

n>0 ai,...,a n ki,...,k n 

x5 kl +...+k„, (10) 



with %Jt ai ...a n {v a — i^ a ; ki, . . . , k n ) being the nth cumulant defined by 

t , ^lnH HS [...] 



^ai...a n ( l/ Q! i^ai ^ 



(11) 



<9<5u; kliai . . . dSu knian 

The nth cumulant 37t ai ... an coincides with the Fourier transform of the n-partical connected 
correlation function of the RS [20]. 5k 1 +...+k„ i s the Kronecker symbol. The nth cumulant 
depends on both the wave vectors ki and the partial chemical potentials u a — iui a . 



Using (fl0l) - (fTTj) we can rewrite ©-((ID as follows 

(66 p) exp | - — Yl ^(k)5p Ka 5p. Kf3 

a,/3 k 

( — l) n 

o.\...ot ri \yct i^Q!) ki : . . . , k n ) 

a k n>2 ai,...,a„ ki,...,k n 

xSu kl>ai . . . Su knta J kl+ ... +kn \, (12) 



where Smf is the GPF of the model in the MF approximation. 

a. Gaussian approximation Now we consider the Gaussian approximation of B.[u a ] set- 
ting dJl ai ...a„ = for n > 3. Then, after integration in (|T2|) over 5u ka we obtain 



a ] - / (d<5p)exp{ - -Y^2^(k)6p k)a 8p^ K py (13) 



o,/3 k 



where C a p(k) is the Fourier transform of the two-particle direct correlation function in the 
Gaussian approximation 



CafiW = + _!=<*£(*). (14) 



Cap(k) is the Fourier transform of the direct correlation function of a two-component hard- 
sphere system. It is connected with Wl a p(k) by the relation 



Cf s (k)Wl 2 (k) = 1, (15) 



where C^ s {k) denotes the matrix of elements C^{k) / ^N a Np and the matrix of elements 
%Raf3(k)- I is the unit matrix. It should be noted that C a p(k) is connected to the ordinary 
direct correlation function c a p(k) by 



C a p(k) = — — c a/3 (k), 

< Pa> 



where p a =< N a > /V. 



n order to determine the CV connected with the order parameter we follow the ideas 
of i3, EJ and introduce independent collective excitations by means of the orthogonal 



transformation 



<5p k ,+ = A(A0£k,i + C{k)£ k;2 

5p k ,_ = B{k)^ k!l +D{k)^ i2 . (16) 



The explicit form of coefficients A(k), B(k), C(k) and D(k) are given in Appendix A. As a 
result, ffT3l) is rewritten as 

E G [v a ] = S MF [z7 Q ~ iw a ] S' A (df) exp | - - y^ea(fe)6caC-k,a}, (17) 

a=l,2 k 

where eigenvalues £i(A;) and ^(fc) are found to be 

e lf2 (*) = \ (c ++ (k)+C-4k) ± [(C ++ (fc) -C__(*)) 2 + ACl_(k)] 1/2 ^j . (18) 

Here we are interested in the gas-liquid critical point. Thus, we are now in a position to 
study equations (fl6j) and ([TBI in the long- wavelength limit. In this case equations (fT6l) have 
the form 

which in turn leads to the relations 

V 1 + z 

£0,2 = 7== (<W + ^Po,-) • (20) 
V l + 2 

Introducing CVs p ,iv = #Po,+ + <Vo,- an d Po,q — z 5po,+ — 5po- that describe long- wavelength 
fluctuations of the total number density and charge density, respectively, we can rewrite (120]) 
in the form 

p _ 1 

^ = 7rT^lTT7^ + TTlH- (21) 

As is seen, CV £o,i describes fluctuations of the charge density. In the general case z ^ 1, 
£0,2 is a linear combination of CVs po.iv and po,Q with the z-dependent coefficients. At z — 1, 
CV £0,2 describes solely fluctuations of the total number density. Thus, we suggest that CV 
£0,2 is connected with the order parameter of the gas-liquid critical point. 
At k = one finds that 

£i(/c = 0) = 00, (22) 



e 2 (k = 0) ' 



(-s^TTi) + 5 ++ (0) + 2V ^" (0) + ^ (0) ) ' (23) 



where T* and 5 are given by (jHJ)-© and p* = pa± is a reduced total number density. 
Equation fl22|) leads to 

G QQ (k = 0)=0, (24) 

where Gqq(/c = 0) is the Fourier transformation of the charge-charge connected correlation 
function; equation ( 1241) reflects the fact that the first moment Stillinger-Lovett rule is satisfied 
in the Gaussian approximation. 

At 5 = 0, e^ik = 0) reduces to the form 

(1 + z) 2 1 

where S^O) is the two-particle structure factor of a one-component hard-sphere system at 
k = 0. In the Percus-Yevick (PY) approximation 29] 

s ^> - (25) 

7T 

where r] = —pa 3 is the packing fraction. It is worth noting that §2(6 = 0; k = 0) takes only 
6 

positive values. It suggests that in the size-symmetric case no phase separation between two 
uniform phases can be found at the Gaussian level of the description. 

Equation £2(6 7^ 0; k = 0) = leads to the gas-liquid spinodal curve in the Gaussian 
approximation 



4irp*8 2 



T* s = (c^(0)/, + 2c^(0)/Vi+c^(0))- 1 . (26) 

Equations (Tl9l) - (l2"6T) are analogous to those obtained in [l9(] but for another type of the 
regularization of the Coulomb potential inside the hard core. The trends of the critical 
parameters calculated from the maximum point of spinodal fT26|) are consistent with the 
corresponding trends found in [19]: at the fixed z the critical temperature T* is a convex 
down function of 5 while the critical density p*(S) is a convex up in 5; both T* and p* increase 
at a given 5 > and decrease at a given 5 < when z increases. Therefore, only some of 
the trends are correctly predicted within the framework of this approximation. In order 
to properly describe the effects of size and charge asymmetry on the critical parameters 
one should take into account the terms of the higher-order than the second order in the 
functional Hamiltonian (Ti~2l . We consider this task below. 
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III. CRITICAL PARAMETERS OF ASYMMETRIC PMS: BEYOND THE GAUS- 
SIAN APPROXIMATION 



In order to study the gas-liquid critical points of asymmetric PMs we use the method 
proposed in [21]. First we pass from the initial chemical potentials u + and v_ to their linear 
combinations 

ZU + -V- V+ + ZV- 

Vl = ^TTW' V2 = 7TT? • (2?) 

Chemical potentials V\ and v 2 are conjugate to CVs £0,1 and £0,2, respectively. Since we 
suggest that CV £0,2 is connected with the order parameter, v 2 appears to be of special 
interest in our study. 



Following the ideas of 
mation (EED 



2l| we start with the logarithm of GPF in the Gaussian approxi- 



ln E G [p a \ = In S HS [v a \ - ^ In det [1 + $ c aJt 2 ] , (28) 



k 



where §c and 9^2 are matrices of elements /30^g(/c) and Wl a p(is a — iuj a ; k), respectively. 

We approximate cumulants Wl a/ 3(k) by their values in the long- wavelength limit putting 
yR a p{k) = Wlap^k = 0) = VJlap. If it is remembered that lnSns and Tl ai a 2 -a n are functions 
of the full chemical potentials we can present z/i and v 2 as 

v x = v\ + X°Au u v 2 = v\ + A°Az/ 2 , 

with v\ and v 2 being the MF values of v\ and u 2 , respectively and Av\ and Av 2 being the 
solutions of the equations 

dlnE G (^ 2 ) =Q; (29) 



= X ((N + ) HS + z(N.) HS ). (30) 



91nS G (i/i,z/ 2; „ 



dAu 2 

We self-consistently solve equations (I29l) -( l30l) for the relevant chemical potential Az/oby 



means of successive approximations keeping terms of a certain order in parameter A 21]. 
To this end, we expand (128]) in powers of Aui and Av 2 

hxEafr, v 2 ) = E °n 1 2j A^r-A^, (31) 

n>0 i„>0 

where 

9 n lnS G (z/i,z/ 2 ) 



Mt\vlA) 



dAv^dAvl 
10 



Ai/i=0,A!/ 2 =0 



The first nontrivial approximation corresponding to Aui = yields the following expression 
for Az/2 

AU2 = 2V W ++ +^2^1 + z 2 DJl ] S det [1 TQdm] (^++ 

+(34>-(k)S 2 + 2(34> + „(k)S 3 ) , (32) 



where 

s 1 = m +++ + zm ++ _, s 2 = 9H + __ + s 3 = mt ++ _ + ,rfm + __. (33) 

Apart from 93t aia2 , formulas ([32l) -( l33l) include the third order cumulants %Jl ai a 2 a 3 or equiv- 
alently the third order connected correlation functions of the RS. 
Finally, we can write the full chemical potential v 2 in the form 

V2 = + 4 + A "2, (34) 

where 

»? S = + (35) 
vl + z 

with z/^ 5 (u^ s ) being the hard-sphere chemical potential of the ath species and being 
the combination of the self-energy parts of chemical potentials u + and z/_ 

4 = ~ 2V J TT ^ E (^?+(*) + */^--(*)) • (36) 
Now some comments are in order: 

• Here we consider the two-component hard-sphere system as a RS. In this case the 
analytical expressions for second order cumulants 9Jt aiQ , 2 can be obtained in the PY 
approximation using the Lebowitz' solution 3l|, [32J. The corresponding formulas for 

(k — 0) are given in Appendix B. 

• In order to derive the expressions for the third order cumulants one can use the recur- 
rent relation 

art -art (0 ) - dma ^- a «-^°> • • j (o 7 \ 

where is the MF value of chemical potential v ai which due to the electroneutrality 
condition coincides with the hard-sphere chemical potential of the aijth species. 

Formulas (132|) - (!36|) will be used for the study of the gas-liquid phase equilibria in asym- 
metric PMs. Below we consider some particular cases. 
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A. Monovalent PMs with size asymmetry 



First we consider a monovalent PM with size asymmetry corresponding to z = 1 and 
A 7^ 1. Because of symmetry with respect to the exchange of + and — ions, only A < 1 (or 
A > 1) need be considered in this case. 

We put z = 1 in (!32|) - (l36|) and consider the PY approximation for the thermodynamic 
and structural functions of the two-component hard sphere system 3l|, |32J. For cumulants 
9JT Ql a 2 and Tl ai0 , 20 , 3 we use formulas (l43]) - (!58l) from Appendix B. The Fourier transforms of 
the interaction potentials are given by (j3j)-([7|). 



The explicit expressions for v^ s and v\ are obtained using the results of Ref. 32J sup- 
plemented by the electroneutrality condition. They are given in Appendix C. 

Based on the expressions fl32l - fl36l) (at z — 1) supplemented by the Maxwell construction 
we calculate the coexistence curves and the corresponding critical parameters for different 
values of A. Estimates of the critical point values of T* and p* are given by their values for 
which the maxima and minima of the van der Waals loops coalesce. The estimated values 
of the critical parameters are presented in table 1. 

TABLE I: Critical parameters T* = ksTa±/ q 2 and p* = p c o\ of the monovalent PM for different 
values of A 



A 


/TT* 
1 C 


10 2 p* 


1.0 


0.0848 


0.907 


0.75 


0.0831 


0.816 


0.5 


0.0786 


0.637 


0.25 


0.0709 


0.433 


0.1 


0.0586 


0.195 



Figures 1 and 2 demonstrate the effects of size asymmetry on the critical parameters 
of the monovalent PM. In Fig. 1 the critical temperature T* depending on A is shown by 
the solid circles for A ranging from 0.1 to 1. As is seen, a qualitative agreement with the 
simulation data shown by the open circles is obtained. The results obtained in the MSA 
are shown by the open squares. In Fig 2 the dependence of the critical density p* on A is 
shown. Similar to the computer simulation findings our results indicate a decrease of the 
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0.04 - 
0.03 o- 




0.2 0.4 0.6 0.8 1.0 

X 



FIG. 1: Critical temperature T* of the monovalent PM as a function of size asymmetry. Open 
circles correspond to the results of simulations 9j; open squares are MSA results via the energy 
route 15[ and solid circles correspond to the results of the CV based theory. 



critical density with the increase of A but the figures obtained in this approximation turn 
out to be too small. 




FIG. 2: Critical density p* c of the monovalent PM as a function of size asymmetry. The meaning 
of the symbols is the same as in Fig. 1. 

It should be noted that the RPM limit turns out to be a special case. When z — A = 1 
expression (I3"2|) reduces to the form that corresponds to the random phase approximation 



(RPA) |2l| 
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B. Equisize PMs with charge asymmetry 



Another particular case is an equisize PM with charge asymmetry corresponding to A = 1 
and z ^ 1. At A = 1 (5 = 0) the expression 



(fc)Si + P4>—(k)S 2 + 2[34> + _{k)S 3 



entering (132]) reduces to the form (see Eqs. (JSD"©, (!32|) - (l33|) and the formulas in Ap- 
pendix B): 



T*x 3 



{(l-zf-(l + zi)S 2 (0)) 



(38) 



where the dimensionless temperature T* is given by Eq. (JH]) under condition a± = a. S 2 (0) 
is given by (125]) . 

As a result, Eqs. (1321 - (136]) read as 

l + z h / 1 (l-z 3 )(l-z) 



Au 2 



,HS 



y/TTz 1 v 

a+z) 



z(l 



,HS 



VTT 



z) 2 + (i + z 2 ys 2 {o) 

l+z 1 



where i\ under conditions (12]) is reduced to the form 



1 f 00 x 2 sin(x) 
T* J or + k* sin(i) 



(39) 
(40) 

(41) 



k* = kd<j, is the Debye number. 

In the PY approximation z/ H5 is as follows 

z 



v 



HS 



\n(rj) + ln(l — rj) + 



\n(z) - ln(l + z) + 



7/(14 - 13t/ + 5t/ 2 ) 



l + v_/ 2il //)•"• 

Taking into account only the first term in (|39|) we arrive at the expression for the chemical 
potential v 2 in the RPA. As is seen from Eq. (T4T]) . Az/ 2 does not include the factor of charge 
asymmetry explicitly in this case. The second term is the correction to the RPA resulting 
from the consideration of the higher-order correlation effects, namely the third order. It 
should be noted that Eq. (1391) is obtained in the approximation which is different from 
that considered in 21| (our notation Az/2 corresponds to Az/jy in 21|). In particular, the 



corresponding formula in 



2l| (Eq. (31) in 21J) includes the third and forth order cumulants. 



Putting z = 1 in (!39]) - fj4T|) we arrive at the chemical potential of the RPM in the RPA 

1 



l>2 



V2 



.HS 



2T* 



71 



(42) 
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which reflects the fact of a special symmetry of this model. In order to go beyond the RPA, 
the higher-order correlations should be taken into account as it was done in [24| . 

Based on (|22D-(|1ID we calculate the coexistence curves and the corresponding critical 
parameters for different values of z. The values of the critical parameters T* and p* for 
different z are shown in Table 2. 

TABLE II: Critical parameters of the equisize PM with charge asymmetry for different values of A 



z 


1 c 




1 


0.0848 


0.907 


2 


0.0640 


0.720 


3 


0.0469 


0.549 




FIG. 3: Dependence of the critical temperature on charge asymmetry for equisize PMs. Solid 
symbols correspond to the results of the CV based theory: circles are the results based on ([39]- (HTJ) , 



triang 
z = l 



24]. Open circles are the results of simulations: 



es are the results from [2]L ■ is for RPM 

l^,z = 2-30],z = 4 6]. The dotted line is the result from the RPA. 



Figures 3 and 4 show trends of T*(z 
simulation data {(], Q. The results of 



and p*{z) obtained from fl39l) - fj4TJ) . together with 



2l| are added for comparison. It should 



}e noted 



that the Carnahan-Starling approximation for the hard-sphere system was used in [21] . The 
critical parameters of the RPM obtained within the framework of this theory but in the 
higher-order approximation are shown by the solid squares (T* = 0.0503, p* = 0.042) [24]. 
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FIG. 4: Dependence of the critical density p* = pa 3 on charge asymmetry. The meaning of the 
symbols is the same as in Fig. 4. 

As in |2l|], the trend of T*(z) obtained from (!39|) - (l4T|) qualitatively agrees with simulation 
findings. Moreover, now the numerical values of T c are much closer to the simulation data 
than those found previously. On the other hand, the trend of p*(z) found from (l39l) - (|4T|) 
is inconsistent with the simulations. It follows from [21] that the correlation effects of the 
higher order than the third order should be taken into account in order to get the correct 
trend of the critical density. 



C. PMs with size and charge asymmetry 

Now we use formulas (I32l) - fl36l) for the study of the size- and charge- asymmetric PMs 
with z = 2 and A ^ 1. As before, in order to calculate the coexistence curves and the 
corresponding critical parameters at different values of A we apply the Maxwell construction. 
Our results for the critical parameters T* and p* are given in Table 3. As is seen, both the 
critical temperature and the critical density decrease with the increase of size asymmetry. 

The dependence of the critical parameters on the size asymmetry is shown graphically in 
Figs. 5 and 6, respectively, along with the results of simulations [h]]. In general, the trends of 
T* and p* with 5 are consistent with the simulation findings: 2:1 systems exhibit a maximum 
in both the critical temperature and the critical density when plotted as a function of size 
asymmetry. Similar to simulations, our results (especially the critical temperature) reveal a 
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pronounced sensitivity to S. However, both the critical temperature and the critical density 
found from fl32l) - fl36l) demonstrate maxima at 5 = (A = 1) while the corresponding maxima 
obtained by simulations are shifted towards nonzero values of S (8 > ). Interestingly, p*(5) 
demonstrates the general shape similar to that obtained for the dumbbell system. As before, 
our values of the critical density are more than an order of magnitude lower than those found 
in the simulations 10j. 



TABLE III: Critical parameters of the (2 : 1) PM for different values of 5 (5 = (A - 1)/(A + 1)) 



8 


■* C 




-0.67 


0.054 


0.319 


-0.5 


0.0587 


0.458 


-0.33 


0.0614 


0.545 


-0.2 


0.0630 


0.619 





0.0640 


0.720 


0.2 


0.0611 


0.616 


0.33 


0.0583 


0.528 


0.5 


0.0553 


0.461 


0.67 


0.0529 


0.436 




FIG. 5: Critical temperature of the (2:1) PM as a function of size asymmetry. Solid circles 
correspond to the results of the CV based theory; open circles are the results of simulations [kJ. 
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0,002 



0,0 0,2 
8 



FIG. 6: Critical density of the (2:1) PM as a function of size asymmetry. Solid circles are the 
results of the CV based theory. Open symbols correspond to the results of simulations [10(] : circles 
are spherical electrolytes; squares are dumbbell electrolytes. 

IV. CONCLUSIONS 

In this paper we have studied the effects of size and charge asymmetry on the gas-liquid 
critical parameters of two-component PMs using the CV based theory. The theory allows 
one to take into consideration the effects of higher-order correlations. On the other hand, 
the well-known approximations for the free energy, in particular DHLL and ORPA, can be 
reproduced within the framework of this theory. Recently, this approach has been used for 
the study of the effects of charge asymmetry on the critical parameters of equisize charge- 
asymmetric PMs. It allowed us to calculate, without additional assumptions (such as the 
presence of the dipoles or the higher-order clusters, for example), the trends of both the 
critical temperature and the critical density with charge asymmetry that qualitatively agree 



with Monte Carlo simulation results 



2l|. 



First, we have studied the Gaussian approximation of the functional Hamiltonian of the 
asymmetric PM. The stability analysis has led us to the trends for the critical parameters 
consistent with those obtained in [lj]]. As was shown in (l9|. only some of the effects 
of size and charge asymmetry can be reproduced correctly at this level of consideration. 
Then, we study the critical parameters of size- and charge-asymmetric PMs taking into 
account the higher-order correlation effects. Following the procedure described in [2l|] we 
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have derived in the first nontrivial approximation an explicit expression for the chemical 
potential conjugate to the order parameter. It includes the third-order correlation functions 
of the reference system. Then, this expression was used to study the three versions of PM: 
an equisize PM with charge asymmetry; a monovalent PM with size asymmetry; a size- and 
(2:1) charge-asymmetric PM. 

In conclusion, within the framework of the same approximation we have obtained the 
trends of the critical temperature and the critical density with size asymmetry that qual- 
itatively agree with the Monte Carlo simulation findings: both T* and p* decrease with 
increasing size asymmetry at the fixed z. As regards the charge asymmetry, the present 
approximation for the relevant chemical potential yields a correct trend of the critical tem- 
perature with z and improves the numerical values of T* when compared with our previous 
findings. Unlike the results of [2]]], this approximation leads to the opposite trend of the 
critical density with charge asymmetry although the variation of p* with z is small. We ex- 
pect that the consideration of the correlation effects of higher order will enable us to correct 
the trends of the critical density with charge asymmetry as well as to improve the numerical 
values of the critical parameters. This issue will be considered elsewhere. 



V. APPENDICES 



A. Explicit expressions for the coefficients A(k), B(k), C(k) and D(k) 



1 + ft 



2' 



B(k) 
D(k) 



a 1 



ft 2 



where 



A(k) 
C{k) 



C—(k) - C ++ (k) ± \ {C ++ (k) - C—{k)f + 4C+_(*) a 



2 ■ 

«2 



«1,2 



2C + _(k) 



B. Expressions for 9H aiQ2 (0) and 9JT Q , lQ , 2a3 (0, 0) in the PY approximation 

Using the Lebowitz' solution of the PY equation [31] for a hard sphere system one can 
obtain the explicit expressions for the Fourier transforms of Ornstein-Zernike (OZ) direct 
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correlation functions c£f (k). In the long- wavelength limit they are of the form 



P+ c^(0) = -27/ + (4a + + 3/3 + + 2 7+ ), 



p_d HS (0) = -2?7_(4a_ + 3/?_ + 2 7+ A~ 



V^pT^(O) = -^{A + J B[10/3+_(45 + 3) + 

+6 7+ _(5a + 4) + 4 7+ (6a + 5)]}, 



where the following notations are introduced 

X + fj\ 3 



V = V+ + V~, V+ 



X-7] 



X- + x + X 3 ' 



X- + x + X 3 ' 



1 A yj X + X 

a = —^i — , h — rj 



2A 



x- = — 

X- + x + X 3 ' p ' 



a+ 



1 



{1 - r] 3 + (t/ + + AV)(^ 2 + 4(1 + 77)) 



-3tj_(l - A) 2 [(l + 77+ + A(l + 77_))(1 - 77 + 37/+) 
+r/ + (l-r/)]}, 

' {A 3 (l - r/ 3 ) + (r/+ + AV)(t/ 2 + 4(1 + rj)) 



A 3 (l -r/) 4 

-3r/ + (l - A) 2 [(l + r/ + + A(l + r/_))(l - 77 + 3r/_) 



+Ar/_(l-r/)]} 



0- 



-6 



^++ + ^(1 + A) 2 A(7 + . 



*7-<£- + 4»7+A- 3 (l + A) 2 ^ + _ 



-6 

-3A(1 - X)(X~ 2 rj + g ++ + T]-g—)g+ 



7+ 

7- 
£4-4- 

9— 

9+- 



\ (v+a+ + A 3 r/_a-) 



7+ 
A 3 



74 



274 



1 - A 



(1 


-v) 2 




1 


(1 


-v) 2 




1 



(l-»7) s 



1 + | + ^_(A-1) 



\ , 37/(1 -A) , 

1 + 77; . — (77+ - 77_) 



4(1 + A) 
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A 



V+V- 



The expressions for S a p(k) = 9Jl a p(k)/ ^N a Np can be found from the OZ equations 
S++(k) 

S-(k) 

s + 4k) 



(l-p + c*S(k))(l-p-cZl{k))- 


P+P- 




l-p + c™(k) 






(1 - p + 5g(k))(l - p-c™(k)) - 


P+P- 


{^(k)) 2 


^ P+ p.c^(k) 







(1 - p+c™(k))(l - P-C HS (k)) - P + P-(c^(k)) 2 ' 
Eqs (T45 T) -( l56l) should be supplemented by the electroneutrality condition. 

The explicit expressions for £ +++ (0,0) = 9Jt +++ (0, 0)/iV + and S ++ _(0,0) 

9Jt ++ _(0, 0)/N + can be obtained from the relations 

'dS++(0)\ 



(53) 

(54) 

(55) 
(56) 



S +++ (0,0) = S ++ (0) 



drj 
dS ++ (0) 



'/+ 



5 ++ _(0,0) 




+7?_S__(0) 



dS ++ (0) 



dr\- 



drj_ 

dS ++ (0) 
drj+ 

v+ 



(57) 



(58) 



The expressions for 971 (0,0) and 9Jt H (0,0) can be obtained replacing indices "+" by 

indices "— " and vice verse. 

Final formulas should be supplemented by the electroneutrality condition. 



C. Explicit expressions for v§ and vJ? s 



We obtain for 



V6 



{1+Xjjz + X) 
4T*AVl + z 2 ' 



In the PY approximation vf s has the form 



MS 



l + z 



Inrj + ln(l — rj) + 



Z i / \ , / T](l + T] + T] 2 ) 

ln(z) - ln(z + A ) + ■ 



3r/ 



1 + Z V ' y "" ' ' (1 - T]) 3 

[2z V (l - A) 2 ((1 + X)(z + A 3 ) + V X(z + A 2 )) 



2(1 + z)(l - t?) 3 (z + A 3 ) 2 
-A(z + X)(z + X 2 )(z + A 3 )(l - V ) 2 - 377(2 + A 2 ) 3 (l - V )]] . 
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